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Abstract 

We incorporate explicit Nystrom methods into the RKQ algorithm 
for stepwise global error control in numerical solutions of initial- value 
problems. The initial-value problem is transformed into an explicitly 
second-order problem, so as to be suitable for Nystrom integration. 
The Nystrom methods used are fourth-order, fifth-order and 10th- 
order. Two examples demonstrate the effectiveness of the algorithm. 

1 Introduction 

In two previous papers we have considered the KKrvQz algorithm for step- 
wise control of the global error in the numerical solution of an initial-value 
problem (IVP), using Runge-Kutta methods P[2]. In the current paper, the 
third in the series, we focus our attention on the use of Nystrom methods in 
this error control algorithm for n-dimensional problems of the form 

y"{x) = {{x,y) (1) 
y{xo) = yo 
y'{xo) = y[). 

Note that f is not dependent on y'. We designate this Nystrom-based al- 
gorithm KKNrvQz, and we will show in a later section how any first-order 
IVP can be written in the form ([1]), so that KKNrvQz is, in fact, gener- 
ally applicable. The motivation for considering this modification to KKrvQz 
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is twofold: most physical systems are described by second-order differential 
equations, and Nystrom methods applied to ([T]) tend to be more efficient 
than their Runge-Kutta counterparts. 

2 Relevant Concepts, Terminology and No- 
tation 

Here we describe concepts, terminology and notation relevant to our work. 
Note that boldface quantities are nxl vectors, except for q:[, I„, FJ^, F^, and 
gy, which are n x n matrices. 

2.1 Nystrom Methods 

The most general definition of a Nystrom method (sometimes known as 
Runge-Kutta- Nystrom (RKN)) for solving ([1]) is 

kp = f + Cphi, Wi + Cphiw[ + /if ^ apgk^ p = 1, 2, m 

m 

Wj+i =Wi + hiW- + hf Yl bpkp = Wi + hiF {xi, Wj) (2) 

p=i 

m ^ 

w'i^i = w'i + hiJ2 bpK- 

p=i 

The coefficients Cp, Upg, bp and bp are unique to the given method. If Upg = 
for all p ^ q, then the method is said to be explicit; otherwise, it is known as 
an implicit RKN method. We will focus our attention on explicit methods. 
In the second line of ([2]), we have implicitly defined the function F. We 
treat w- as an 'internal parameter'; for our purposes here, we do not identify 
w' with y', because f is not dependent on y'. The symbol w is used here 
and throughout to indicate the approximate numerical solution, whereas the 
symbol y will be used to denote the exact solution. We will denote an RKN 
method of order r as RKNr and, for such a method, we write 

w[^i=w[ + /i,F^(x„w[,wr). (3) 
The stepsize hi is given by 

hi = Xi_^\ Xi 

and carries the subscript because it may vary from step to step. It is known 
that RKNr has a local error of order r + 1 and a global error of order r, just 
like its Runge-Kutta counterpart RKr. 
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2.2 IVPs in the form y" = f (x, y) 

Consider the n-dimensional IVP 



This gives 



y'{x) = g{x,y) (4) 
y(2;o) = yo- 



// % i^^ y) dyi ^ dgj (x, y) 



where yi is the ith component of y, and Qi is the ith component of g. Clearly, 
we have 

for all j = 1,2, . . . ,n, and so we can write 

y" (x) = f {x, y) . 
The initial values for this second-order problem are then given by 

y{xo) = yo 

y'(a;o) = g(a;o,yo) = yo- 

Hence, any first-order IVP can be transformed into an IVP of the form ([1]). 
This is ideally suited to the Nystrom methods, which are specifically designed 
for this type of IVP. They are also more efficient than their Runge-Kutta 
counterparts; for example, the methods to be used later, RKN4 and RKN5, 
require three and four stage evaluations, respectively, as opposed to RK4 and 
RK5, which require at least four and six stage evaluations, respectively. 

2.3 Error Propagation in RKN 

It can be shown [3] that, for RKr, 

A[+i = w[+i - yi+i = £[+1 + a[A[ (5) 
a[ ^ In + KFlix,,^,), (6) 

where e^^-^^ = O is the local error, A[_,_j^ is the global error and is 

the Jacobian (with respect to y) of the function F*" (xj, w[) associated with 
RKr. The term hiFy {xi, in the matrix cx.^ arises from a first-order Taylor 
expansion of F*^ (xj, Wj) = F*" (xj, y, + A^) with respect to yi. 
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For a Nystrom method RKNr, we have F'^ = F*" (xj, w[) and so, as above, 

Oil = I„ + hiF'y {xiXi) 1 

where is an appropriate constant. Hence, the global error in RKNr is also 
given by ([5]). 

2.4 RKrvQz 

We will not discuss RKrt'Qz in detail here; the reader is referred to our pre- 
vious work where the algorithm has been discussed extensively. It suffices 
to say that KKrvQz uses RKr and KKv to control local error via local ex- 
trapolation, while simultaneously using KKz to keep track of the global error 
in the RKr solution. Such global error arises due to the propagation of the 
KKv global error. KKrvQz is designed to estimate the various components 
of the global error in RKr and RKt; at each node and, when the global er- 
ror is deemed too large, a quenching procedure is carried out. This simply 
involves replacing the RKr and KKv solutions with the much more accurate 
KKz solution, whenever necessary, so that the RKr and RKt; global errors 
do not accumulate beyond a desired tolerance. 

2.5 RKNri;Qz 

The algorithm KKNrvQz is nothing more than RKrvQz with RKr, KKv 
and RKz replaced with RKNr, RKNf and RKNz. Of course, KKl^rvQz is 
applied to problems of the form ([1]), whereas KKrvQz is applied to problems 
of the form (j4]). 

We also report on a refinement to the algorithm: in RKrf Q2;, if the global 
error at is too large, we replace w[ with w^" and then recompute w[_^j^ and 
wJYi, using as input for both RKr and RK^;. This is the essence of the 
quenching procedure. However, in retrospect it seems quite acceptable to 
simply replace w[^^ and wj^i with wf^^; this avoids the need for recomputing 
w[^^ and wJYi, which improves efficiency and, after all, it is the global error 
in w[^_]^ and wj^i, not w[ and w^, that is too large. Both approaches are 
effective, although one is more efficient than the other. It is the more efficient 
approach that we have employed in KKNrvQz. 

3 Numerical Examples 

It is not our intention to compare methods or algorithms but, for the sake 
of consistency, we will apply KKNrvQz to the same examples that we con- 
sidered in our previous work on KKrvQz. In our calculations, we use RKN4, 
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RKN5 and RKNIO which gives the algorithm RKN45Q10. RKN4 and RKN5 
are taken from Hairer et al [S], and RKNIO is from Dormand et al [4J. 
The first of these is the scalar problem 



y' 
?/(0) 



In 1000 



100 / 



which transforms to 



y" 

2/(0) 
2/'(0) 



In 1000 Y 
100 ) 
1 

In 1000 
100 



Solving this problem with RKN45 and RKN45Q10 with a tolerance of 10~^° 
on the absolute local and global errors gives the error curves shown in Figure 
1. The global error obtained with RKN45 is clearly larger than the desired 
tolerance on most of the interval, despite local error control via local extrap- 
olation. However, RKN45Q10 yields a solution with a global error always 
less than the tolerance - the maximum global error in this case is 9.1 x 10~^^. 
The points on the x-axis where this global error decreases sharply correspond 
to the quenches carried out using RKNIO. 

The second example is the simple harmonic oscillator 



y'l 
y'2- 



y(o) 



2/2 

-2/1 


1000 



which has solution 



2/1 {x) 
2/2 {x) 



1000 sinx 
1000 cos a; 



and becomes, in explicit second-order form. 



y(o) 



y'i 

y'i " 



1000 



-2/1 
. ~2/2 . 

. y' (0) 



f (2^,y) 

■ 1000 " 
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Since the solution oscillates between —1000 and 1000, there are regions where 
the solution has magnitude less than unity - here, we implement absolute 
error control - and regions where the solution has magnitude greater than 
unity, where we implement relative error control. With an imposed tolerance 
of 10~* on the local and global errors (relative and absolute) we found a 
maximum global error of ~ 4 x 10^^ in each component when using RKN45, 
and a global error no greater than 0.99 x 10^* with RKN45Q10, on a; G 
[0, 200] . A total of 20 quenches were needed. 

4 Conclusion 

We have considered the use of Nystrom methods in KKrvQz, wherein a com- 
bination of local extrapolation and quenching result in stepwise global error 
control in numerical solutions of IVPs. Two examples have demonstrated 
the success of RKN45Q10. 
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Figure 1: Global errors for RKN45 and RKN45Q10 appUed to the 
exponential test problem. 



